en 



PLh 



c^ 



COMPUTING A LOGARITHM OF A UNITARY MATRIX WITH 

GENERAL SPECTRUM 

TERRY A. LORING 



O, 

^N , Abstract. We analyze an algorithm for computing a skew-Hermitian logarithm of a 

f^ ' unitary matrix. This algorithm is very easy to implement using standard software and 

^ . it works well even for unitary matrices with no spectral conditions assumed. Certain 

examples, with many eigenvalues near —1, lead to very non-Hermitian output for other 
basic methods of calculating matrix logarithms. Altering the output of these algorithms 
to force an Hermitian output creates accuracy issues which are avoided in the considered 
algorithm. 

.^ ' A modification is introduced to deal properly with the J-skew symmetric unitary 

h^ , matrices. Applications to numerical studies of topological insulators in two symmetry 

^^ ' classes are discussed. 

a 
P.: 

1. Introduction 

"^ . While all invertible matrices have a logarithm, indeed many logarithms, the unitary 

,__! ! matrices have the nicest logarithms. Every unitary matrix U has a skew-Hermitian 

lO ' logarithm. Indeed, when U is unitary, the conditions 

\0'. -TV < K <'K, e'^ = U 

^ I specify K uniquely. Still working abstractly, we may take advantage of the finiteness 

^ ■ of the spectrum of U when discussing additional symmetries. For each U there is a 

complex polynomial so that K = p{U) which means, for example, that when U is J- 
skew-symmetric then so must be K. 
k> ' We may regard U* = U~^ as a symmetry, but the non-linearity of the unit circle 

J^ ■ inevitably makes this symmetry a little different from that of being Hermitian {A = A*) 

or symmetric {A^ = ^)- In the case of a matrix being Hermitian, if numerical errors 
start to creep in, we simple replace A by \A+ |yl* in 0{p?) time and have again A* = A 
exactly. The unitary part of the polar decomposition can be used when U*U is only 
close to / but it takes more that 0{n^) time to compute and the result is not exactly 
unitary. For this reason we need to consider errors in ||?7*f/ — /|| a bit larger than machine 
precision. 

Remark. We are using U* to denote conjugate-transpose of a matrix. This would be 
denoted with a dagger in physics. For just the conjugate we use U. 
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What is called self-dual [9] in physics is called J -skew- symmetric [3] or skew Ham,ilton- 
ian p] in computer science. 

Finally, what applied mathematicians call matrix functions are called, in pure math, 
applications of the functional calculus. 

The algorithm discussed here arose in a numerical study in condensed matter physics 
[lOj. In that situation ||f/*?7 — /|| ~ | was typical. Such a matrix U is still very well 
conditioned, so one would expect almost any algorithm for computing a logarithm to 
perform well, at least most of the time. However the study in question had a time- 
reversal symmetry that, through a variation on Kramers pairs, caused the approximate 
unitaries in question to have multiplicity at least two in every eigenvalue. Otherwise the 
approximate unitaries were free to have arbitrary spectrum within the unit circle. Some 
reflection makes one realize that whatever branch of logarithm one uses, a degenerate 
eigenvalue on the branch point has the potential to cause trouble. This is true whether 
or not the degeneracy in the spectrum is caused by a symmetry or by very bad luck. 

We are not claiming it is hard to find K with e'^ ~ U. We are interested in achieving 
e*^ ^ U and K* = K aX the same time. Moreover we want an algorithm that can be 
easily analyzed and also easily modified to take into account additional symmetries. 

To be consistent, we consider the branch of logarithm that has the imaginary part of 
log(A) in the interval (— tt, vr], so in particular log(— 1) = m. We are mainly concerned 
with the operator norm 

\\M\2 



= sup 

VT^O ||V||2 

We also will use the Frobenius norm, i.e. the un-normalized Hilbert-Schmidt norm 



\A\\f 



\ 



j,k=i 



2 

ik\ ■ 



Recall the bounds \\A\\ < \\A\\p < y/n\\A\\. 

One easy solution, that often works well for logarithms of well-conditioned matrices, 
starts by diagonalizing U via an invertible. 

Algorithm 1. 

(1) Compute a diagonalization, W invertible and T diagonal with WTW~^ ~ U. 

(2) Create a diagonal unitary matrix D via Djj = Tjj/ \Tjj\. 

(3) Compute i^o = W\og{D)W-\ 

(4) Output: H = \H* + \Hq. 

This fails badly in special cases, even for small matrices. Given 

-1 
0-1 
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W 



D 



we might consider the approximate diagonahzation 

-1 

exp(— 7ri + Oi) 

to be rather good if 6* > is small, even if a > is rather large. Indeed 

— 1 a(l + exp(— 772 + 6*^)) 



1 a 
1 



WDW~ 



U 



exp(— Tri + 6i) 

a 



U 



which has norm 



However, if we set 



= (1 + exp(-7r2 + ei)) 

|1 + exp(— 7ri + 9i)\ yl + a^ 
Ko = -iW \og{D)W-^ 
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Tx —Tia + 9a 
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-vr 


+ e 







] 







TX 





-2'Ka + 9a 

-71 + 9 



we note that this is not close to the actual logarithm, which is diagonal. It is not close to 
any branch of logarithm applied to U. When we enforce the symmetry hj K = ^Kq + ^Kq, 
we obtain 

71 —7ia + ^9a 

-71 + 9 



K 



-7xa + ^9a 



When a = 1 we have 
and yet 

where 

Working numerically, we find 



\im\\WDW'^-U\ 



71 



-7xa 



—na 



y -71 



Ki 



7T 



lim lle'^ - U\ 

9^0 " ' 



-71 



JK-i 



9 



-7T 
-7T 







^Ki 



-U 



1.2114 



which is dramatically off. 

Examining the methods for computing matrix exponentials, both good and bad, in fiE\ , 
we notice methods based on the Schur decomposition. In theory, the Schur decomposition 
of a normal matrix will result in a an upper triangular factor that is actually diagonal. 
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So almost normal should lead to an almost diagonal factor, but this is too naive an 
approach. However we are not concerned with general almost normal matrices just now, 
only the nice special case of almost unitary matrices. A naive approach using the Schur 
decomposition will work well in this specialized situation. 

The logm function in MATLAB is optimized to so that H = logm(f/) leads to expm(if ) 
being very close to U, even when U is badly conditioned. It is not optimized to convert 
the approximate relation U*U ^ I into the approximate relation H* ?» —H. The Schur- 
Parlett algorithm [6J behind logm and f unm is primarily intended for entire functions. It 
can behaive unexpectedly when eigenvalues are at and near branch points. 

Algorithm 2. 

(1) Compute logarithm Hq of U via logm. 

(2) Output: H = \H* + \Hq. 

A far more ambitious project would be to find algorithms for matrix functions that, 
given almost normal matrices as input, lead to almost normal matrices on output. The 
subject of almost normal matrices is discussed in a recent survey [5] and is full of subtle 
problems. 

2. SCHUR FACTORIZATION OF NEAR UNITARY MATRICES 

We define the deviation from unitary to be the number ||[/*f/ — /||. 
In finite precision arithmetic, we can expect the deviation from unitary to almost never 
equal zero. Not surprisingly this error is easy to handle. 

Lemma 1. Suppose U is in M„(C) and Uv = Av for some unit vector v and scalar A. 
Then 

||A|^-l| < \\U*U-I\\. 

Proof. Let P denote the orthogonal projection on the one-dimensional space spanned by 
A. Then 

||A|^-l| = ||(|A|^-l)P|| 

= \\P{U*U -I)P\\ 

< \\U*U-I\\. 

D 

Remark 2. Notice then 

\^T^^W*u~^\ < \M < \^i + \\u*u -i\\ 

so that as long as \\U*U — /|| < 1 we have 

||A|-1| < \\U*U -I\\. 

Lemma 3. Suppose T in M„(C) is upper triangular and let D denote the diagonal matrix 
corresponding to the diagonal ofT. Then 

\\T -D\L < J2(n-1) \\T*T-Ip . 

II 1 1 i^ V \ / 1 1 II 
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Proof. We know that the Tjj are all eigenvalues of T so 






1 < iiT*r-/|| . 



Every element on the diagonal of a matrix has absolute value at most the norm of that 
matrix. Applied to T*T — I this tells us 



\T*T-I\\ > 



-i + Ei^' 



kj\ 



k=l 



i-1 

^ \ ^ IT' l'^ I IT 1^ II 

k=l 



'JJ\ 



SO 



J-1 



^iTfcjf < 2||T*r-/||. 



fc=i 



Summing these error bounds we learn 



\T - DWpY < 2{n - 1) \\T*T - I\\ . 



D 



Remark 4. We get here the estimate 



IIT-Dllp < ^/2{n - 1) \\T*T - I 
which we can compare to Henrici's estimate [11] 



I^-^IIf< 



n — n 
12 



Irri^rTi T~'T~'*II 2 



As expected, we are getting a better estimate for almost unitary matrices than works in 
the more general case of almost normal matrices. 

Lemma 5. Suppose T in M„(C) is upper triangular with || T*T — /|| < 1. If we let D 
denote the diagonal matrix with 

1 



D 



n 



T ■ 

IT I ■'■^ 



^JJ\ 



then 



\\T - D\\ < ( V2(^ - 1) + 1 ) 11^*^ - ^11 ' • 
Proof. The previous lemma shows that with E diagonal and Ejj = Tjj we have 



\T - E\\ < \\T - E||p < ^/2(n - 1) \\T*T - IP 



Dealing with diagonal matrices is easy, and we find 

1 



\D — E\\ = sup 



T, 



jj 



T.. 

IT. .!"'•' 



'J3\ 



< \\T*T - I\ 
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by Remark m SO 



\T - D\\ < ^2{n-l) \\T*T - /|p + \\T*T - I\ 



< ( V2(r2-l) + lj ||T*T-/||5. 

D 

Theorem 6. Suppose U is in M„(C) and that unitary Q and upper-triangular matrix T 
are a Schur factorization for U, meaning U = QTQ* . If we define D to be the diagonal 
unitary matrix with 

D = T 

I ^ jj I 



\U-QDQ*\\ < U2{n-l) + l] \\U*U - Ip . 



\T*T - I\\ = \\Q* {U*U - I) Q\\ = \\U*U - I\ 



then 

Proof. We have 

and 

\\T - D\\ = \\Q* {U - QDQ*)Q\\ = \\U - QDQ*\\ 

so this foUows immediately. D 

3. DiAGONALIZING MATRICES THAT ARE CLOSE TO UNITARY 

We may face a matrix U where U*U is as close to / as can be expected within the 
constraints of machine precision. For Hermitian matrices, the eigensolvers in LAPACK 
produce eigenvectors that are "always nearly orthogonal to working precision" lU §4.7.1]. 
No comparable promise is made in other eigensolvers, but off-the-shelf algorithm such as 
ZGEES in LAPACK computing a Schur factorization form a good substitute. We get a 
simple algorithm for finding a unitary eigensolver for a unitary matrix. 

Algorithm 3. 

(1) Compute a Schur factorization, Q unitary and T upper-triangular, QTQ* ^ U. 

(2) Create a unitary diagonal matrix D via Djj = Tjj/ \Tjj\. 

(3) Compute Qlog(L')Q*. 

If ||f/*f/ — /|| is a larger than machine precision, then we can compute the unitary part 
of U and proceed as above. For simplicity of programming, we will compute the polar 
decomposition via the singular value decomposition. We are not advocating this method, 
heeding the warnings in p!2]. 



Algorithm 4. 

(1) Set V = U{U*U)-k 

(2) Compute a Schur factorization, Q unitary and T upper-triangular, QTQ* ^ V. 

(3) Create a unitary diagonal matrix D via Djj = Tjj/ \Tjj\. 

(4) Compute giog(L')C*; 
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deviation 
from unitary 



Backwards Error lie ^^ — U 



n 



\U*U -I\ 



Algorithm 1 



Algorithm 2 



Algorithm 3 



Algorithm 4 



Algorithm 5 



4.11082e-15 



0.12285 



0.12340 



4.30902e-15 



4.41189e-15 



4.13976e-15 



16 



5.02961e-15 



0.04407 



0.04465 



6.31606e-15 



6.47906e-15 



6.13171e-15 



32 



6.33082e-15 



0.09286 



0.09099 



9.36391e-15 



9.30067e-15 



8.99073e-15 



64 



1.10432e-14 



0.01952 



0.01641 



1.38378e-14 



1.39124e-14 



1.32675e-14 



128 



1.34734e-14 



0.01999 



0.02239 



2.29980e-14 



2.32533e-14 



2.26790e-14 



256 3.19324e-14 



0.06131 



0.06158 



4.49885e-14 4.31729e-14 4.42639e-14 



Table 1 . We compare algorithms that produce matrices H with H* = H 
exactly but with e^^ ~ U with varying accuracy. These are run on unitary 
matrices with at least two eigenvalues near —1. For each n the data shown 
are averaged over 30 test matrices. The noise variable in the code in the 



Appendix was set to 10" 



-15 



n 



-0.56 



The exponent —0.56 was determined to 



keep the deviation from unitary close to constant across different matrix 
sizes. 



Our focus is on U with deviation from unitary in the range (0, 0.1). Newton's method of 
approximating the unitary part of U is very effective in this situation. Newton's method 
here sets V = U and iterates the replacement V ^ ^ (V + (y~^)*). For our purposes, it 
makes sense to use two interations. 

Algorithm 5. 

(1) Setl^i = i(f/+(f/-i)*). 

(2) SetV=^{V, + {V,-'y). 

(3) Compute a Schur factorization, Q unitary and T upper-triangular, QTQ* ^ V. 

(4) Create a unitary diagonal matrix D via Djj = Tjj/ \Tjj\. 

(5) Compute Qlog(D)Q*. 

One could compute when to stop the iterations for best accuracy, following, fT2]. As 
simpler methods work here, where we have such well-conditioned matrices, we do not 
pursue this option. 

Even one iteration of Newton's method has advantages. 

Lemma 7. Suppose U is in M„(C) and \\U*U — I\\ < |. // 



then 



and 



v=-,{u+{u-r) 



\V*V-I\\ < \\U*U -I\ 



\U-V\\ < \\U*U -I\ 
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deviation 
from unitary 



Time 



n 



\U*U-I\ 



Algorithm 1 



Algorithm 2 



Algorithm 3 



Algorithm 4 



Algorithm 5 



4.11082e-15 



0.00013s 



0.00190s 



0.00009s 



0.00014s 



0.00016s 



16 



5.02961e-15 



0.00034s 



0.00618s 



0.00024s 



0.00036s 



0.00037s 



32 



6.33082e-15 



0.00139s 



0.02104s 



0.00104s 



0.00146s 



0.00145s 



64 



1.10432e-14 



0.00877s 



0.05853s 



0.00719s 



0.00926s 



0.00886s 



128 



1.34734e-14 



0.06111s 



0.13750s 



0.05442s 



0.06843s 



0.06218s 



256 3.19324e-14 



0.33397s 



0.73867s 



0.29425s 



0.35710s 



0.32437s 



Table 2. Timing data for the algorithms and data set from Table[TJ Times 
were computed in a non-parallel environment. 



Proof. Let U = DQ be the alternate-side polar decomposition of U, so Q is unitary and 
D is positive semi-definite. Then 



|f/*f/-/|| = 11^*1)2^-/11 = ||L)2_j|| = gup |a2-1| 

AG(7(D) 



and 



and 



\V*V -I\ 



\u-v\\ 



1{D + D-^) 



D-l{D + D-^) 



sup 

xecT{D) 



i{x+x-r-i 



sup |A- i (A + A~^)| . 

AG(7(D) 



This lemma reduces to routine algebra, showing that 



1 , V7 

2 - - 2 



implies 



and 



< A^-l 



i(A + A-) -1 

A-|(A + A'i)|<|A2-l| 



D 



Theorem 8. Suppose U is m M„(C) with \\U*U — I\\ < | and let 



V = ^ ([/+([/-)•). 
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deviation 
from unitary 



Backwards Error lie ^^ — U 



n 



\U*U ~I\ 



Algorithm 1 



Algorithm 2 



Algorithm 3 



Algorithm 4 



Algorithm 5 



1.19608e-05 



0.33294 



0.33294 



9.42901e-06 



5.98042e-06 



5.98042e-06 



16 



1.27157e-05 



0.18364 



0.18364 



1.01912e-05 



6.35784e-06 



6.35784e-06 



32 



1.25555e-05 



0.24113 



0.24113 



1.01824e-05 



6.27775e-06 



6.27775e-06 



64 



1.21902e-05 



0.21304 



0.21304 



9.96141e-06 



6.09511e-06 



6.09511e-06 



128 



1.18934e-05 



0.24809 



0.24809 



9.74603e-06 



5.94671e-06 



5.94671e-06 



256 



1.15296e-05 



0.25532 



0.25532 



9.45829e-06 5.76481e-06 



5.76481e-06 



Table 3. As in Table [T], but with the noise variable set to 10 ^n '^•^^. 



// there is a unitary Q and an upper-triangular matrix T so that V = QTQ* , and if we 
define D to be the diagonal unitary matrix with 

D = T 

I -^ jj I 

then 

\\U -QDQ*\\ < U2{n~l) + 2) \\U*U -I\\. 

Proof. We have 

\\V-QDQ*\\ < U2{n-l) + l) \\V*V - I\\^ < U2{n - 1) + l) \\U*U - I\\ 

and 

\\u -V\\ < \\U*U -I\\. 

D 

It is worth keeping in mind that the spectral decomposition of the polar part of U 
leads to the theoretical best unitary diagonalization, with D diagonal, Q unitary and 



l|f/-g*Dg|| =max||A-i| xeaUu*uy- 

and the best general lower bound is 

\\U-Q*DQ\\ >^\\U*U-I\\. 

Lemma 9. Suppose U is m M„(C) and \\U*U - /|| < |. // 

1 



and 



I4 = 5(c+(f-')-) 
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deviation 
from unitary 



Backwards Error lie ^^ — U\ 



n 



\U*U -I\ 



Algorithm 1 



Algorithm 2 



Algorithm 3 



Algorithm 4 



Algorithm 5 



3.98986e-01 



0.48072 



0.48083 



2.91776e-01 



1.87829e-01 



1.87829e-01 



16 



3.99794e-01 



0.76608 



0.76629 



3.00304e-01 



1.88836e-01 



1.88836e-01 



32 



4.01659e-01 



0.74454 



0.74577 



3.02891e-01 



1.86026e-01 



1.86026e-01 



64 



4.04556e-01 



0.81147 



0.81217 



2.99941e-01 



1.85677e-01 



1.85677e-01 



128 



3.92340e-01 



0.86934 



0.87020 



2.90461e-01 



1.79970e-01 



1.79970e-01 



256 3.81114e-01 



1.23535 



1.23575 



2.84104e-01 1.75207e-01 1.75207e-01 



Table 4. As in Tabled! but with the noise variable set to 0.3n ^■^^. 



then 
and 



\V*V-I\\ < ^\\u*u 



\u-v\\ < ^\\U*U-I\\. 



Proof. Tracking the spectrum of positive parts as before, we are looking at 



1 A^ + l 

2 A 



+ 1 A^ + 6A2 + 1 



' (1^) 4A(A2 + 1) 

and all that is needed is some algebra and calculus to verify that 



2 - - 2 



implies 



and 



A^ + 6A2 + 1 
4A (A2 + 1) 

A^ + 6A2 + 1 



< 



4A(A2 



1^ 



25 



- 10 ' 



1 . 



D 



As before, we use this to get an estimate on the algorithm that uses two iterations of 
Newton's method followed by a Schur decomposition. 



Theorem 10. Suppose U is in M^ 



with \\U*U - /|| < I and n>3. Let 



v^=-,iu+{u-r) 
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and 



v=l{v^+{v,-r) 



If there is a unitary Q and an upper-triangular matrix T so that V = QTQ* , and if we 
define D to be the diagonal unitary matrix with 



D 



n 



IT I ■'■^ 



■-Jjl 



then 



Proof. We have 



and 



\U-QDQ*\\ < ^y^\\U*U - if + ^ \\U*U - I\ 



10 



\V*V-I\\ < —\\U*U-lf 
I 11-25 " " 



\V - QDQ*\\ < ( x/2in - 1) + 1 ) \\V*V - Ip 



i^/2{n- 


-1+1) 


\u*u- 


-I 




lovnii 


U*U- 


-i\'y 









and 



\u -v\\ < ^\\U*U -I\ 



u 



Rather extreme input data is need to highhght the advantage of Algorithm [51 We test 
algorithms [2ll5] on unitaries, and approximate unitaries, that have multiple eigenvalues 
very near —1. See Tables [T]|ll 

Not much can be said regarding the error in the output H of Algorithm [5] vs. the "true 
logarithm." The trouble is that the set of "exactly connect" answers jumps around in a 
most discontinious fashion. If a, /3 are in (— tt, tt] with a ~ — vr and /3 ?^ vr then for any 
unitary Q, consider the unitary 



U = Q 







e'^ 
The only exactly correct output given U on input is 



Q*. 



Q 



a 

/3 



Q* 



and this greatly depends on Q. 
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deviation 
from unitary 


Backwards Error e *^ 


-U 


n 


\\U*U-I\ 


Algorithm lA 


Algorithm 2A 


Algorithm 6 


8 


2.96904e-15 


1.17188 


1.12520 


3.27683e-15 


16 


3.27269e-15 


1.10402 


1.11826 


4.50363e-15 


32 


4.12604e-15 


1.13947 


1.15103 


6.68904e-15 


64 


8.33263e-15 


0.77344 


0.77214 


1.00208e-14 


128 


1.02179e-14 


1.46559 


1.46102 


1.52540e-14 


256 


2.46375e-14 


1.19984 


1.19987 


2.78177e-14 



Table 5. We compare algorithms that produce matrices H with H* = 
H'^ = H exactly but with e*^ ~ U with varying accuracy. These are run 
on unitary, J-skew-symmetric matrices with at least four eigenvalues near 
— 1. For each n the results shown are averaged iver 30 test matrices. The 



noise variable was set to 10 



-15 



n 



-0.56 



in the code in the appendix. 



4. J-SKEW-SYMMETRIC UNITARIES 

When dealing with Hamiltonians for systems with certain time-reversal symmetry, we 
need an extra involution on matrices, the dual. Working in N-hj-N blocks, the dual 
operation on n-hj-n matrices with n = 2N is defined as 



A B 
C D 






A" 



It is easy to check this obeys the same axiom as the transpose. In particular it commutes 
with the adjoint. See [16], for example. The dual operation is not unique, but is fixed 
once we specify 

/ 
-/ 



J 



Then X^ = —JX^J. A matrix X is J-skew-symmetric when X^ = X and so if and only 
if (XJ)^ = -XJ. 

There is a J-skew-symmetric Schur decom,position result for J-skew-symmetric matri- 
ces, a.k.a. the skew-Hamiltonian Schur decomposition result from [2]. 



Theorem 11. Let n = 2N. If X'^ = X in M„(C) then there is a unitary Q with Q'^ 
and a matrix S with X = QSQ* and so that 



Q* 



S 



T B 

TT 



where T is upper triangular. Moreover there is an 0{n^) algorithm to find Q and T. 
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deviation 
from unitary 



Time 



n 



\U*U -I\ 



Algorithm lA 



Algorithm 2A 



Algorithm 6 



2.96904e-15 



0.00019s 



0.00223s 



0.00250s 



16 



3.27269e-15 



0.00034s 



0.00503s 



0.00616s 



32 



4.12604e-15 



0.00124s 



0.01752s 



0.01492s 



64 



8.33263e-15 



0.00773s 



0.05169s 



0.03858s 



128 



1.02179e-14 



0.05879s 



0.16380s 



0.14007s 



256 



2.46375e-14 



0.34309s 



0.52146s 



0.79723s 



Table 6. Timing data for the algorithms and data set from Table [5 



Proof. This is the complex analog of [21 §2.1], discussed in detail in [TOl §9.1]. One uses 
some variation on the Paige / Van Loan algorithm [21 §2.1], involving careful combinations 
of Givens rotations and partial Householder reflections to create a unitary Qi with Q\ = 
Ql and so that X = QiSiQ\ with 

TT 
WTW* now finishes the job, as 

W _ 

W 



Si 



An ordinary Schur decomposition Ti 

Q 

has the needed symmetry, while 

S = 
has the correct block structure. 



Qi 



w 



w 



Si 



w 



w 



n 



For comparision purposes, we say we are using Algorithm [2]A, etc, if we add a final 
step to the algorithms above that replaces H by \HK 

Notice that S becomes upper-triangular if we just reverse the order of half the basis 
elements. 

Theorem 12. Let n = 2N. Suppose U is a J -skew-symmetric matrix in M„(C) and 
that unitary Q and matrix T form a J -skew-symmetric Schur decomposition for U. If 
we define D to he the diagonal unitary matrix with 

1 



D 



n 



T.. 

IT ^ " 



^n\ 



then D is J -skew- symmetric, diagonal and 

\\U-QDQ*\\ < fv^2(r2-l) + l) \\U*U - I\ 
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deviation 
from unitary 


Backwards Error e *^ 


-U 


n 


\U*U-I\ 


Algorithm lA 


Algorithm 2A 


Algorithm 6 


8 


6.80865e-06 


0.50219 


0.50219 


3.40432e-06 


16 


7.91995e-06 


0.27931 


0.27931 


3.95997e-06 


32 


8.21901e-06 


0.32397 


0.32397 


4.10950e-06 


64 


8.30081e-06 


0.43814 


0.43814 


4.15041e-06 


128 


8.15713e-06 


0.30607 


0.30607 


4.07856e-06 


256 


8.02983e-06 


0.38732 


0.38732 


4.01491e-06 


ABLE 


7. As in Tabled but with the noise variable set to 10~^n~^'^^ 



If we approximate the polar part of U by Newton's method, the symmetry f/" = [/ is 
preserved at each iteration, since 

{I {u + (£/-)•))' = lu' + (([/-')•)• = lu* + (([/-')')• = it/' + (([/')-)• . 

Theorem 13. Let n = 2N. Suppose = JJ is m M„(C) with \\U*U -I\\<j and let 



and 



v,=\{u+(u-r) 



v = ^(v. + (vr')') 



// a unitary Q and matrix T form a J -skew- symmetric Schur decomposition for V , and 
if we define D the diagonal unitary matrix with 

1 



D 



n 



m 



:T. 



331 



33 \ 



then D is J -skew-symmetric and 



\\U-QDQ*\\ < J^^\\U*U - if + ^\\U*U - I\\ . 

Theorem [13] gives theoretical justification for that the following algorithm produces 
approximately correct output. This is very similar to the algorithm used in [17]. Since 
this is built out of known algorithms, we don't have anything to say it is 0{n^). 

Algorithm 6. 

(1) SetV, = l{U+iU~r). 

(2) Set V='^{V^ + {¥,-')*). 

(3) Compute a J- skew- symmetric Schur factorization, as in Theorem [T3l with QTQ* ^ 
V. 

(4) Create a unitary diagonal matrix D via Djj = Tjj/ \Tjj\. 
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deviation 
from unitary 


Backwards Error e *^ — f/ 


n 


\U*U-I\\ 


Algorithm lA 


Algorithm 2A 


Algorithm 6 


8 


2.12731e-01 


0.45404 


0.45282 


1.05752e-01 


16 


2.43013e-01 


0.55252 


0.55205 


1.18737e-01 


32 


2.51907e-01 


0.83382 


0.83331 


1.21165e-01 


64 


2.61957e-01 


0.96758 


0.96743 


1.24065e-01 


128 


2.62333e-01 


0.94305 


0.94325 


1.23593e-01 


256 


2.56967e-01 


1.03853 


1.03868 


1.21172e-01 



Table 8. As in Tabled but with the noise variable set to 0.3n' 



-0.56 



(5) Compute Qlog(L')Q*. 

For data on timing and accuracy, see Tables [SHHl These tables were produced with 
the MATLAB code listed in the appendix. That code has almost no optimization. In 
particular, it does not take advantage of the symmetries Q" = Q and D'^ = D that hold 
at the top of the loop in the Paige - Van Loan algorithm. 



5. Logarithms in Physics 

In quantum mechanics, it is standard to exponential a skew-Hermitian operator to get 
a unitary, as this is how one moves from the Hamiltonian to the time evolution operator. 
When then is a periodic time-dependent Hamiltonian Ht of period T, the definition of 
quasi-energy depends on the Floquet Hamiltonian Hp defined via 

It does not matter if the principal branch of logarithm is used to define Hp but it 
is important in numerical studies that Hp be computed to be Hermitian. The study of 
Floquet topological insulators jTl] is an important special case of a system with a periodic 
time-dependent Hamiltonian. Some Floquet topological insulators can only be explained 
by keeping track of a form of time- reversal symmetry |l3i equation 21]. In that case one 
has a self-dual Floquet Hamiltonian. 

The logarithms of unitary matrices arise in another way in physics, in particular in 
the study of the more typical topological insulators where there is a time-independents 
Hamiltonian. For finite lattice models of non-interacting fermions, the dimension and 
symmetry class determine if distinct topological phases can occur. Physically these phases 
are ordinary insulators and topological insulators [H IT9] . 

Joint working with Hastings [10] established that the ordinary insulating phases can be 
characterized by the existence of localized vectors that form a basis of low-energy space 
(localized Wannier functions) where the basis preserves an appropriate symmetry. We 
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explained in Section 4 of [10] how to translate this question into a question about almost 
normal, almost unitary, or almost commuting matrices. 

The All symmetry class is one where there is a certain flavor of time-reversal invariance. 
What this means mathematically is that one starts with J-skew-symmetric and Hermitian 
matrices for the Hamiltonian and the position observables. There is a promising method 
for computing the spin Chern numbers that involves self-dual logarithms. This was 
introduced in [17]. The formula in ii'-theory used there is now validated by the theorems 
in \15\ and the results in [15] . 
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7. Appendix 

7.1. Code for General Unitaries. 

function [unitary_error, time, accuracy] = testLogs(n, noise,NumbNewt) 
time = 0; 
accuracy = 0; 

/ / / y / / / y o / o / o / o / o / o / o / o / o / o / o / o / o / o / o y o / o / o / o / o y oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% Create a test matrix, approximately unitary 

% Not fully random as we want some spectrum near -1 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% exponentiate a random self-adjoint matrix to create a unitary 

K = 0.25*( rand(n)+i*rand(n)-rand(n)-i*rand(n)) ; 

K = K + K' ; 

K = (4*pi /norm(K))*K; 

"/oexponentiate i*K to get a unitary 

Q = expm(i*K) ; 

% create a "random" diagonal, to form the test unitary 
D = diag(exp(2*pi*i*[0.5, 0.5, rand(l,n-2)] )) ; 
U = Q*D*Q' ; 

Zadd noise 

U = U+noise*(rand(n)+i*rand(n) -rand(n)-i*rand(n)) ; 

unitary_error = norm(U'*U - eye(size(U))) ; 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% Get a logarithm by diagonalizing U = W*D*inv(A) , then 
% force it to be hermitian 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

tstart = tic; 

[W,D] = eig(U); 

X = diag(imag(log(diag(D)))) ; 

H = W*X*inv(W) ; 

H = (1/2)*(H + H'); 

time(l) = toc(tstart); 
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accuracy(l) = norm(expni(i*H)-U) ; 

/ / / y / y o / o y o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o y oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o In /o /o /o /o In la la la la la la la la la la la la la la 

"la Get a logarithm using logm, then force it 
"la to be hermitian 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la 

tstart = tic; 

H = (-i)*logm(U); 

H = (1/2)*(H+H'); 

time (2) = toe (tstart); 
accuracy(2) = norm(expm(i*H)-U) ; 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la 

"la Get a logarithm by Schur decomposition then 
"la force it to be hermitian 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la 

tstart = tic; 

"la compute the Schur factorization 

[q,T] = schur(U); 

D = diag(imag(log(diag(T)))) ; 

H = Q*D*Q' ; 

H = (1/2)*(H + H'); 

time (3) = toe (tstart); 
accuracyO) = norm(expm(i*H)-U) ; 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la 

"la Get a logarithm from the "exact' ' unitary part 
"la followed by Schur decompositionthen 
"la force it to be hermitian 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la la 

tstart = tic; 

[A,Z,B] = svd(U); 

V = A*B' ; 

"la compute the Schur factorization 

[Q,T] = schur(V); 

D = diag(imag(log(diag(T)))) ; 

H = Q*D*Q' ; 

H = (1/2)*(H + H'); 

time (4) = toe (tstart); 
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accuracy(4) = norm(expm(i*H)-U) ; 

/ y / y o / o y o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o / o y oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% Get a logarithm by from an approximate unitary part 
% followed by Schur decompositionthen 
% force it to be hermit ian 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

tstart = tic; 

V = U; 

for round = l:NumbNewt 

V = 0.5*(V + inv(V) '); 
end 

% compute the Schur factorization 

[Q,T] = schur (V); 

D = diag(imag(log(diag(T)))) ; 

H = Q*D*Q' ; 

H = (1/2)*(H + H'); 

time (5) = toe (tstart); 
accuracy(5) = norm(expm(i*H)-U) ; 

end 

7.2. Code for J-skew symmetric Unitaries. 

function [unitary_error, time, accuracy] =testLogsDual(n, noise,NumbNewt) 

% force n to be even 
n = n + mod(n,2) ; 
N = n/2; 

time = 0; 
accuracy = 0; 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% Create a test matrix, approximately unitary 

% Not fully random as we want some spectrum near -1 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o / /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o / /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o / /o /o 

% exponentiate a random anti-self -adjoint and self -dual matrix to 

% create a symplectic unitary 

K = 0.25*( rand(2*N)+i*rand(2*N)-rand(2*N)-i*rand(2*N)) ; 

K(N+1:2*N,N+1:2*N) = -K(l :N, 1 :N) . ' ; 

K = (1/2) *(K - dual(K)); 
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K = (1/2)*(K + K'); 

K = (4*pi /norm(K))*K; 

Zexponentiate i*K to get a symplectic unitary 

Q = expm(i*K) ; 

% create a "random" diagonal, for form the test unitary 
D = exp(2*pi*i*[0.5, 0.5,rand(l,N-2)] ) ; 
D = diag([D,D]); 
U = q*D*Q' ; 

Zadd noise, but keep self -dual 

U = U+noise*(rand(2*N)+i*rand(2*N)-rand(2*N)-i*rand(2*N)) ; 

U = (1/2)*(U + dual(U)); 

unitary_error = norm(U'*U - eye(size(U))) ; 

/ / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / y / / / 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% Get a logarithm by diagonalizing U = W*D*inv(A) , then 
% force it to be hermit ian and self dual 

/ y / / / y o / o / o / o y oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

tstart = tic; 

[W,D] = eig(U); 

X = diag(imag(log(diag(D)))) ; 

H = W*X*inv(W) ; 

H = (1/2)*(H + H'); 

H = (1/2) *(H + dual(H)); 

time(l) = toc(tstart); 
accuracy(l) = norm(expm(i*H)-U) ; 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

% Get a logarithm using the default logm, then force it 
y„ to be hermitian and self -dual 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

tstart = tic; 

H = (-i)*logm(U); 

H = (1/2)*(H+H'); 

H = (1/2)*(H + dual(H)); 

time (2) = toe (tstart); 
accuracy(2) = norm(expm(i*H)-U) ; 

oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 
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% Get a logarithm using Newton's method 

% followed by a stuctured Schur decomposition, then force it 

% to be hermitian and self -dual 

/ y / y o / o / o / o y o / o / o / o y oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy oy 
/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o In /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o In /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 

tstart = tic; 

y„Move toward a unitary with Newton's method 

V = U; 

for round = l:NumbNewt 

V = 0.5*(V + inv(V) '); 
end 

% compute the Schur factorization 

% First the Paige - Van Loan algorithm 

[Q,D] = PVL(V); 

D = D(1:N,1:N); 

% compute the little Schur factorization, finalize Q and D. 

[Q1,T1] = schur (D); 

Q(1:2*N,1:N) = Q(l : 2*N, 1 :N)*Q1 ; 

Q(1:2*N,1+N:2*N) = Q(l :2*N, 1+N: 2*N)*conj (Ql) ; 

D= [diag(imag(log(diag(Tl))))] ; 

D = [D, zeros (N); zeros (N), D.']; 

H = Q*D*Q' ; 

H = (1/2)*(H+H'); 

H = (1/2)*(H + dual(H)); 

time (3) = toe (tstart); 
accuracyO) = norm(expm(i*H)-U) ; 

end 



% The Paige - Van Loan algorithm 
function [Q,D] = PVL(V); 

n = size(V) ; 
n = n(l)/2; 

Q = eye(2*n) ; 

D = V; y„ Expect to keep V = Q*D*Q' 

for k=l: (n-1) 

"/oHouseholder to fix most of the bottom half in column k 
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if k<(n-l) 

[v, beta] = gallery Chouse ' ,D(k+l+n:2*n,k)) ; 
Zon left of D, top, then bottom 
D(k+l:n,l:2*n) = D(k+l:n,l:2*n) ... 

- (conj (beta)*conj (v))*(conj (v) ' *D(k+l :n, 1 :2*n)) ; 
D(k+l+n:2*n,l:2*n) = D(k+l+n:2*n,l:2*n) ... 

- (beta*v) * (v ' *D (k+l+n : 2*n , 1 : 2*n) ) ; 
"/oon right of D, top, then bottom 
D(l:2*n,k+l:n) = D(l :2*n,k+l :n) ... 

- (beta*(D(l:2*n,k+l:n)*conj (v) )) *conj (v) ' ; 
D(l:2*n,k+l+n:2*n) = D(l:2*n, k+l+n :2*n) ... 

- (conj (beta)*(D(l:2*n,k+l+n:2*n)*v)) *v'; 
Zon right of Q 

Q(l:2*n,k+l:n) = Q(l:2*n,k+l:n) ... 

- (beta*(Q(l:2*n,k+l:n)*conj (v) )) *conj (v) ' ; 
Q(l:2*n,k+l+n:2*n) = Q(l:2*n, k+l+n :2*n) ... 

- (conj (beta)*(Q(l:2*n,k+l+n:2*n)*v)) *v'; 
end 

ZA symplectic Givens rotation to clear out D(k+l+n,k) 

[G,y] = planerot([D(k+l,k) ;D (k+l+n, k)]) ; 

%on left of D 

top.row = D(k+l,l:2*n) ; 

bottom_row = D (k+l+n, 1 :2*n) ; 

D(k+l,l:2*n) = 0(1 , l)*top_row + G(l,2)*bottom_row; 

D(k+l+n,l:2*n) = G(2, l)*top_row + G(2,2)*bottom_row; 

Zon right of D 

left_col = D(l:2*n,k+1) ; 

right.col = D(l :2*n, k+l+n) ; 

D(l:2*n,k+l)=conj (G(l, l))*left_col+conj (G(l,2))*right_col; 

D ( 1 : 2*n, k+l+n) =conj (G(2 , 1) ) *lef t_col+conj (G(2 , 2) ) *right_col ; 

%on right of Q 

left_col = Q(l:2*n,k+1) ; 

right.col = Q(l :2*n, k+l+n) ; 

Q(l:2*n,k+l)=conj (G(l, l))*left_col+conj (G(l,2))*right_col; 

Q(l:2*n,k+l+n)=conj(G(2,l))*left_col+conj(G(2,2))*right_col; 

ZHouseholder to fix top half in column k 
if k<(n-l) 

[v, beta] = gallery ( 'house ' ,D(k+l :n,k)) ; 
Zon left of D, top, then bottom 
D(k+l:n,l:2*n) = D(k+1 :n, 1 :2*n) ... 
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- (beta*v)*(v'*D(k+l:n,l:2*n)) ; 
D(k+l+n:2*n,l:2*n) = D(k+l+n:2*n,l:2*n) ... 

- (conj (beta)*conj (v))*(conj (v) ' *D(k+l+n:2*n, 1 :2*n)) ; 
Zon right of D, top, then bottom 

D(l:2*n,k+l:n) = D(l :2*n,k+l :n) ... 

- (conj (beta)*(D(l:2*n,k+l:n)*v)) *v' ; 
D(l:2*n,k+l+n:2*n) = D(l:2*n,k+l+n:2*n) ... 

- (beta *(D(l:2*n,k+l+n:2*n)*conj (v)))*conj (v) ' ; 
Zon right of Q 

Q(l:2*n,k+l:n) = Q(l:2*n,k+l:n) ... 

- (conj (beta)*(Q(l:2*n,k+l:n)*v)) *v' ; 
Q(l:2*n,k+l+n:2*n) = Q(l:2*n,k+l+n:2*n) ... 

- (beta*(Q(l:2*n,k+l+n:2*n)*conj (v)))*conj (v) ' ; 
end 

end 

end 



function Y = dual(X) 

N = size(X); 
N = N(l)/2; 

Y(N+1:2*N,N+1:2*N) = X(l :N, 1 :N) . ' ; 
Y(1:N,1:N) = X(N+1 :2*N,N+1 : 2*N) . ' ; 
Y(1:N,N+1:2*N) = - X(l :N,N+1 : 2*N) . ' ; 
Y(N+1:2*N,1:N) = - X(N+1 :2*N, 1 :N) . ' ; 

end 
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